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Abstract — We consider the problem of recovering signals from 
their power spectral density. This is a classical problem referred 
to in literature as the phase retrieval problem, and is of 
paramount importance in many fields of applied sciences. In 
general, additional prior information about the signal is required 
to guarantee unique recovery as the mapping from signals to 
power spectral density is not one-to-one. In this paper, we assume 
that the underlying signals are sparse. 

Recently, semidefinite programming (SDP) based approaches 
were explored by various researchers. Simulations of these 
algorithms strongly suggest that signals upto o(^n) sparsity 
can be recovered by this technique. In this work, we develop 
a tractable algorithm based on re weighted -minimization that 
recovers a sparse signal from its power spectral density for 
significantly higher sparsities, which is unprecedented. We also 
discuss the limitations of the existing SDP algorithms, and 
provide a combinatorial algorithm which requires significantly 
fewer measurements to guarantee recovery. 

Index Terms — Phase Retrieval, Semidefinite Programming 
(SDP), Reweighted /i -minimization 

I. Introduction 

In many practical measurement systems, the power spectral 
density of the signal, i.e. the magnitude square of the Fourier 
transform, is the measurable quantity. Phase information of 
the Fourier transform is completely lost, because of which 
signal recovery is difficult. This problem occurs in many 
areas of engineering and applied physics, including X-ray 
crystallography [[3], astronomical imaging J?], microscopy, 
optics lis], blind channel estimation and so on. 

Recovering a signal from its Fourier transform magnitude, 
or equivalently its autocorrelation, is known as phase re- 
trieval. The mapping from signals to their Fourier transform 
magnitude is not one-to-one, and hence unique recovery is 
not possible in general. Additional measurements or prior 
information about the signal is required in order to uniquely 
recover the underlying signal. Consttaints on the signal's 
values like non-negativity, bounds on the signal's support 
(locations where the value is non-zero), and more recently 
sparsity ||9|, are commonly used as prior information. 
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Considerable amount of research has been done over the 
last few decades (H], lH) and a wide range of heuristics have 
been proposed (see ||6l), a comprehensive survey of which 
can be found in Q. El provides a theoretical framework to 
understand the heuristics, which are in essence an alternating 
projection between a convex set and a non-convex set. Such 
methods often converge to a local minimum, hence drastically 
reducing the chances of successful signal recovery. 

Recently, the phase retaeval problem was recast as a semi- 
definite programming problem (see ifTO) . ifTTI and ifTZI ). In 
[10], additional measurements with different illuminations, 
which is possible in an optical setup, are used to make unique 
recovery feasible. In ifTTI . lfT2l and lITSl . the underlying signals 
are assumed to be sparse, which is a reasonable assumption 
in applications like X-ray crystallography, microscopy and 
astronomical imaging. 

Numerical simulations of the existing techniques based on 
SDP strongly suggest that signals upto o{y/ri) sparsity can be 
recovered with an arbitrarily high probability. This behavior 
for the phase retrieval problem was rigorously explained in 
ifTTl . ifTSl . |fT9l . II20II consider the "generalized" phase retrieval 
problem and observe a similar behavior In this work, we 
develop an algorithm based on the idea of reweighted li 
minimization to solve the phase retrieval problem for signif- 
icantly higher sparsities. We also provide certain theoretical 
guarantees and discuss the limitations of the SDP based tech- 
niques, and develop a combinatorial technique which requires 
far fewer measurements to guarantee recovery. 

The remainder of the paper is organized as follows. In 
Section 2, we formulate the phase retrieval problem and recast 
it as a SDP problem. In Section 3, we discuss the limitations of 
the existing SDP-based techniques and develop an algorithm 
based on reweighted li minimization in Section 4. In Section 
5, we develop a measurement system using a combinatorial 
approach which requires far fewer measurements to guarantee 
recovery. Section 6 presents the results of the numerical 
simulations and concludes the paper. 

II. Problem Formulation 

Suppose X — (zq, xi, ....a;„_i) is a real-valued discrete- 
time signal of length n and sparsity k, where sparsity is defined 
as the number of non-zero entries. Let y = {yo,yi, ■■■yn-i) 



be its Fourier transfonn, i.e., 

y = Fx (1) 

where F is the nxn DFT matrix. The phase retrieval problem 
can be mathematically stated as 

find X 

subject to |y| = |Fx| (2) 

Since magnitude square of Fourier transform and autocor- 
relation are Fourier pairs, the phase retrieval problem can be 
reformulated as recovering signals from their autocorrelation, 
i.e., 

find X (3) 

subject to fli — XjXjj^i < i < n ~ 1 
j 

where a ~ {ao,ai, ....,a„_i) is the autocorrelation of x. 

Observe that the operations of time-shift, flipping and global 
sign-change do not affect the autocorrelation, because of 
which there is a trivial ambiguity. The signals resulting from 
these operations are considered equivalent, and in all the 
applications it is considered good enough if any equivalent 
signal is recovered. 

The problem (|2]l is hard to solve due to non-convex con- 
straints. We can relax the constraints into a set of convex 
constraints by embedding the problem in a higher dimensional 
space, a technique popularly known as lifting. Note that (|2]l 
contains n constraints of the form \yi \ = If/^xj, where f, is the 
i*^ column of F. Squaring both sides, the constraints can be 
rewritten as j^jp = |x-'"f*fj^x|. Suppose we define X = xx-'", 
the problem can be recast in terms of X as 

find X 

subject to |y,p = <race(MiX) 0<i<n-l 

ranfc(X) = 1 & X ^ (4) 

where M, = f*ff . 

In terms of the autocorrelation a, the lifted problem can be 
formulated as 

find X 

n 

subject to ^ ~ o-i < i < n ~ 1 

ranfc(X) = 1 & X ^ (5) 

III. Background 

The program (|4) can be reformulated as a rank minimization 
problem as follows 

minimize ranfc(X) 

subject to |y,p = trace(MiX) < i < n - 1 

X ^ (6) 

This is a non-convex problem as rank is a non-convex function. 
It has been shown in lfT4l that trace minimization is the 



tightest convex relaxation of rank minimization for positive 
semidefinite matrices. This relaxation is not useful in the 
phase retrieval setup as irace(X) corresponds to the energy 
of the signal x, which is fixed by the magnitude of the 
Fourier transform. [fl5l proposes log-determinant function as 
a surrogate for rank in such problems, i.e., 

minimize log det (X + el) (7) 
subject to |y,|^ = trace(M,X) 0<i<n-l (8) 
X^O (9) 

This heuristic tries to minimize a concave function in a convex 
domain, which can be done using gradient descent approach. 
This method was explored for the phase retrieval setup in llTOll . 
ifm . Simulations suggest that the algorithm converges to a 
rank 1 solution with high probability if the underlying signal 
is o{y/n) sparse. 

In |12|, we explored a two-stage recovery process to prov- 
ably solve In the first stage, we use information about 
the support of the autocorrelation to recover the support of 
the signal (see iflTl ). In the second stage, we solve the SDP 
with known support ifTSl . It was empirically observed and 
theoretically shown that signals were recovered with arbitrarily 
high probability if the sparsity was o{^/n). However, if the 
support information was available by other means, it was 
observed that the program recovered signals up to roughly 
o(n) sparsity. 

IV. Recovery Algorithm 

In this section, we develop an iterative algorithm based 
on reweighted -minimization to solve the phase retrieval 
problem outside the o{y/n) sparsity regime. 

A two-stage approach like |l2j wouldn't work as the support 
of the autocorrelation becomes full if the signal has sparsity 
greater than 0{^Jnlog{n)). Trace minimization in the phase 
retrieval setup has two issues: trivial ambiguities have same 
objective function, trace is fixed because of which we will be 
solving a feasibility problem only. Weighted li minimization 
( [Tol l intuitively overcomes these issues and promotes sparse 
solutions. 

minimize trace (V|X|) 

subject to |yjp = trace(MjX) 0<z<n-l 

X)pQ (10) 

where V is a weight matrix, which can be designed to promote 
the necessary structure in the solution. 

Simulations suggest that ( fTOl ) has rank 1 solutions with 
high probability if the sparsity is o{\/n) and V is chosen 
randomly, but fails outside the o{^/n) region. However, the 
largest eigenvalue turns out to be considerably stronger than 
the other eigenvalues, and the eigenvector corresponding to it 
happens to contain a lot of information about the support of 
the signal even significantly outside the o{^/n) region, which 
is not very surprising. 



This strongly suggests the possibility of an iterative algo- 
rithm, which at every iteration also knows "a lot" about where 
the signal's non-zero entries can be. Algorithm 1 uses this 
information by doing a re-weighted minimization at every 
iteration. The weights corresponding to prospective support 
locations are set to zero to encourage the signal to choose 
those locations in the next iteration, and the weights outside 
this region is chosen to be positive. 



Algorithm 1 Phase Retrieval Algorithm 

Input: The magnitude of the Fourier transform |y|, maximum 

number of iterations 

Output: The underlying signal x 

1. Initialize V by by choosing its entries from [0, 1] uniformly 
at random 

2. Solve the optimization problem 

minimize trace(V|X|) 

subject to |yip = trace(MiX) < i < n - I 

X)^0 (11) 

3. If rank{'X.) = 1, return X, else calculate Xi = xixf, 
where Xi is the best rank-1 approximation of X 

4. Update V as follows: Yij if |xi| and |xj| are greater 
than a certain threshold, choose the remaining entries from 
[0, 1] uniformly at random 

5. Iterate until convergence or maximum number of iterations 

6. Calculate Xi and return xi 



V. Phase Retrieval as a compression problem 
Theorem IV. II was proved for the phase retrieval problem in 

CD, im. 

Theorem V.l. Signals can be recovered from their power 
spectral densities up to time-shift, reversal and global sign 
with probability \ — 5 for any S > if 

1) k = o(V^) 

2) k entries are chosen uniformly at random 

3) n> n{S) 

We can instead fix the sparsity k and consider the phase 
retrieval problem with partial power spectral density informa- 
tion (for example 1221 ). In particular, one might have access 
to only certain frequencies. To solve this problem, we will use 
some classical results in the compressed sensing literature. 

Theorem V.l (Candes, Romberg, Tao, El]). Let F be the DFT 
matrix. Consider the random matrix A obtained by choosing 
m rows of F uniformly at random. If n. is a k sparse vector, 
X can be recovered from observations Ax with arbitrarily 
high probability if m, > O(fclogn) via the following ii 
minimization: 

min ll^lli (12) 

X 

subject to Ax = Ax 



The following theorem gives a useful result for the recovery 
of a sparse signal from its partial power spectral density. 

Theorem V.3. Let x be a k sparse vector whose support is 
chosen uniformly at random. Suppose we observe its power 
spectral density at m distinct frequencies chosen uniformly 
at random. If m > O(fc^logn), x can be recovered from its 
partial power spectral density with arbitrarily high probability. 

Proof: If X is a fc-sparse signal, its autocorrelation can 
have at most fc^ non-zero entries. Since power spectral density 
is the Fourier transform of the autocorrelation, using Theorem 
IV. 2 1 whenever m > 0{k^ \ogn,), the autocorrelation of x can 
be recovered from partial power spectral density observations 
via £i minimization with arbitrarily high probability. ■ 
While Theorem IV. 3 1 gives a tractable algorithm for the 
recovery of a from partial power spectral density observations 
s = Aa, one can consider a more sophisticated approach . The 
program (fT2t tries to solve for a fc^-sparse solution without 
using the extra information that the resulting signal should 
be a valid autocorrelation. We propose the following method 
which might be of interest for future directions. 

min ||a||i (13) 

a 

subject to Aa = s 

n 

Xj j-i-i = ai < i < n — 1 

ranfc(X) = 1 & X ^ (14) 

Observe that the only nonconvex constraint in ( fTJt is (fT4l l. 
We believe that solving ( fT3l ) can substantially increase the 
performance and instead of O(fc^logn) measurements, just 
0(fclog7i) measurements might suffice which is similar to 
that of typical compressed recovery of a fc sparse vector. 
However, further investigation is required to relax (UTi in a 
useful manner as opposed to dropping it. 

Overall, Theorem IV. 3 1 is subject to a o(Y^)-sparsity bot- 
tleneck since the full power spectral density corresponds to n 
measurements. While we do not provide theoretical guarantees 
for Algorithm 1, when full power spectral density is available, 
our algorithm seemingly beats the o(Y^)-sparsity bottleneck. 
In section [VT] we see that the recoverable sparsity is much 
higher than o{^/n). 

A. Two-stage recovery 

In lfT2l . we try to solve the sparse phase retrieval problem 
via a two-stage approach where the first step involves finding 
the support of the signal from the support of the autocorrela- 
tion. We will now argue that such an approach is inherently 
subject to the o{y/ri) bottleneck, as there is no way of finding 
the support of the signal when its sparsity is greater than 

Lemma V.l. Suppose x is a k-sparse signal whose support 
is chosen uniformly at random, and whose nonzero entries 
are continuous i.i.d. random variables. Then, there exists a 



constant c such that whenever k > c^nlog{n), support of 
the autocorrelation is full with arbitrarily high probability. 

Proof: Without loss of generality, we can assume that 
each location belongs to the support of the signal with prob- 
ability - independently as the same proof will apply for k- 
sparse signals with standard modifications. 

For a particular distance d, if no two non-zero entries in 
the signal are separated by a distance d, we can say that 
d doesn't belong to the support of the autocorrelation. This 
probability can be bounded by (1 — fc^/n^)"/^ which is upper 
bounded by e-'="/2" Union bound tells us that the probability 
of the support of the autocorrelation not being full is less 
than ne~^ which goes to zero if A; > c^J nlog{n) for 
sufficiently large n. ■ 

B. Relation to Gaussian Phase Retrieval 

Our results on partial power spectral density can be related 
to the "generalized" phase retrieval problem, where the ob- 
servations are of the form |gfx| for i.i.d. complex standard 
normal vectors {gj}i!^i- While this problem is structurally 
similar to phase retrieval, it is considerably simpler as there 
are no trivial ambiguities like time-shift and flipping. 

Assuming x is a sparse vector, ||T9| and ll20l analyze the 
following semidefinite program for the tractable recovery of 
X up to a global phase ambiguity. 

min ||X|U + Al|Xl|i (15) 

X 

subject to {^g,gf , = (gjgf , xx"^) l<i<m 

Naturally, one would wish that the unique minimizer of ([15} 
to be xx^ to be able to recover x up to phase ambiguity. 
(flST l tries to capture both sparsity and low rankness of the 
underlying matrix xx^. Interestingly, both |fT9l and ll20l 
suggest that as long as m < 0(min{fc^, n}), recovery using 
( flSl l is not possible with very high probability for any choice 
of regularizer A. To summarize, even if k^ <^ n, one would 
still need r2(fc^) measurements for recovery, which indicates 
an o{^/n) bottleneck for generic Gaussian measurements too. 

On the other hand, Ull, ||20l and JIT] show that, if one 
is able to search for a low rank and sparse matrix, then, 
xx-^ can be recovered with only k log ^, measurements which 
illustrates a significant performance gap between tractable 
recovery and intractable combinatorial search. Interestingly, 
lfT2l illustrates a similar gap for the phase retrieval problem 
with Fourier measurements. 

Overall, we have seen that for the two important class of 
phaseless measurements (Fourier and Gaussian), the current 
theoretical guarantees for the existing algorithms are subject 
to a strong o(Y^)-sparsity bottleneck. A natural question is 
whether it is possible at all to do sparse phase retrieval in a 
tractable way with small number of measurements. 

C. Combinatorial Approach 

To address this question, we consider the problem of re- 
covering X from phaseless measurements by making use of a 



specific choice of measurements. In particular, we show that 
one can tractably recover a /c-sparse signal from phaseless 
measurements by using only O(fclogn) measurements with 
very high probability. This shows that phase retrieval with 
optimal number of measurements (up to a log n factor) is in 
fact possible. 

Theorem V.4. Suppose x is an arbitrary k-sparse vector 
where k > 2. Let {zi}™ be i.i.d. vectors with i.i.d. entries 
distributed as: 

with probability 1 — 1/ k 

(16) 

A/ (0, 1) with probability 1/k 

Let {a^jbi}™^]^ be i.i.d. vectors with i.i.d. entries distributed 
as exp(i0) where 9 is uniformly distributed in [0, 2-k). Denote 
the function R" x R" — >■ R" that returns entrywise products 
of two vectors by ■. Assume, we observe the measurements: 

|(a, •z„x)|M(b, •z„x)|2 (17) 

for I < i < m. 

Then, xx* can be recovered with high probability in 0(mn) 
time, whenever m > ck log n for some constant c > 0. 

Remark: In general, we don't need the knowledge of 
sparsity k for the design of the measurement operator 
This can be handled by introducing an extra factor 
of logn in the required number of measurements, i.e. 
klog^{n) measurements. The reason is we can take cfclogn 
measurements designed for each of the sparsity levels of 
ki = 2* for I < i < [logTi], i.e. use fe; instead of k in ( fT6b . 
Then, k will lie between ki and ki+i for some i and the 
same proof argument will work. 

The proof of Theorem IV.4I can be found in the appendix. 
While our proposed recovery algorithm requires a perfectly 
sparse signal, we emphasize that, our intention with Theorem 
IV.4I is demonstrating the possibility of tractable recovery 
rather than coming up with a robust algorithm for artificially 
designed measurements (fTSI l. 

VI. Numerical Simulations 

In order to demonstrate the performance of Algorithm 1, 
numerical simulations were performed for different values of 
signal length n and sparsities k. 

For a given n and sparsity k, the k support locations 
were chosen from the n locations uniformly at random. The 
signal values at the support locations were chosen from [0,1] 
uniformly at random. The Fourier transform magnitudes of the 
signal were computed and provided as input to Algorithm 1 . If 
the output signal matched the input signal, it was counted as a 
success, else it was counted as a failure. Two other semidefinite 
program based algorithms (CandesPR: ifTO) . HassibiPR: [il2J ) 
and the traditional alternating projection algorithm (GS: 161) 
were used for comparison purposes. 

The results of the numerical simulations are shown in Figure 
1. The probability of successful recovery is plotted against 
various sparsities for n = 32 and n = 64. 100 simulations 



were performed for each sparsity to get the average success 
probability, which was calculated as the percentage of success- 
ful recovery. The existing SDP based algorithms start fading 
at around o{y/n) sparsity, the figure clearly demonstrates 
the superiority of using reweighted minimization outside the 
o(-y/n) sparsity region. 

n.32 
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Fig. 1 . Success rate of recovery for n = 32 and various sparsities 
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Fig. 2. Success rate of recovery for n = 64 and various sparsities 
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Appendix 



Proof of Theorem [ 

Proof: The proof provides a tractable recovery algorithm 
and consists of two main steps. 
Support recovery: Denote the support of x by S* C [n] 
and the support of by Si and let Cj = • a^. Now, 
consider the inner product c*x. Since nonzero entries of 
have continuous distribution \f S C^ Si ^ %, |c,*x| 7^ almost 
surely. Hence whenever a measurement |c*xp = we can 
deduce that 5 n 5"^ = 0. Let: 

I = {l<i<m\Sir\S (18) 

Clearly, {S'ij's are i.i.d. supports and for each i we have: 



P(S'nS', = 0) = (1-1/A:)' 



1 

> - 

- 4 



(19) 



Following ( fT9l ), by the law of large numbers, w.h.p. |/| > 
m/8. Conditioned on |/| > m/8, the probability that j e S 
is not contained in IJ^gj Si is at most (1 — l/fc)'^' < (1 ^ 

Assuming m > 8fc log n and using a union bound: 



\s % U 



< exp(- — ) < 71-1 (20) 
8k 



which will approach 0. Hence, the exact support of the signal 
can be found by simply taking the union of sets Si satisfying 
c*x = and then complementing it (0{mn) time). Next, 
with the knowledge of support, we proceed with the recovery 
of X up to an overall phase ambiguity. 
Signal recovery: Recovery of the signal given its support will 
be performed in two steps. We first show that magnitudes of 
nonzero entries of x can be found. 
Recovering magnitudes: Assume S'i n S* is a singleton j e 
S. Since we already have the knowledge of S from previous 
part, we can immediately deduce: \xj\^ = ^^|^"^|2 ■ Then, 
we simply need to ensure that w.h.p., for all j G S there 
exists 1 < i < m satisfying S D Si = j. Probability of this 
not happening for a fixed j E S is: 

(l4(l4)-r<exp(-|). (21) 

After union bounding over all j G S, whenever m > 
8fclogn, all I a 
l-fcexp(-t^ _ 
Recovering the relative phases: Next, we consider the 
measurements satisfying \Si D S\ ^ 2, 1 < i < m. For 
fixed i we have: 



j\ can be found with probability at least 

) > 1 -n-i. 



n\s^ n ^1 



(1 



fc-2 



> 



1 



(22) 



Overall, w.h.p. there are m/10 measurements satisfying 
\Si n S"! = 2. Let: 



I = {l<i<m\[j\S^^\S\ = 2} 



(23) 



Next, form the k vertex graph obtained by connecting the 
nodes j, I whenever {j, 1} = Si OS for some i G I {0{mn) 



time). Observe that, each i <E I picks an edge in this 
graph uniformly at random. Overall, we perform at least ^ 
picks with replacement. This graph is connected with high 
probability whenever the chance of an edge being picked 
is more than , ll24ll . This happens when m > cfclogfc 
for c > 10 since: 



P(edge is picked) > 1 - (1 - ■ 
> 1 — exp(- 



\m/10 



5k{k- 1) 



> 1 



clogfc 
exp( ^) 



;logfc 



(24) 



(25) 



(26) 



5fc ' 5k 

Now, w.h.p. the graph is connected. Find a spanning tree 
T in this graph, which can be done in O(k^) time, ||25]| . 
Set the phase of the initial node to 0. Then, the rest of the 
phases of the nodes are uniquely determined as follows. 
Recovering relative phases: To begin, consider an edge of 
T between nodes j, I where {j, 1} = SiD S for some i. Fix 



i.e. phase. From measurements: 



I (aj • z,,x) \ ,\> 
and with the knowledge of 



(27) 

we can find: 



?R.{zjZia*aiXjXi) and ^{zjZib*biXjXi) where 5R(-) returns 
the real part of a number. This information is equivalent to: 
^{ax),^{bx) where a = a*ai, b = b*bi have uniformly 
random phases and I . We will argue that the phase 

of X is uniquely determined when ^{ax),3i{bx),a,b,\x\ 
are known. Assume a = exp(i0i), b = exp(i02) and 
X — \x\ exp(i6'). Then, 



^{ax) = |a;|cos(6' + 6li) 
^{bx) = |x|cos(6' + 6l2) 



(28) 
(29) 



gives us two linearly independent equations (almost surely) 
and two unknowns (sin 0, cos 6*). Overall, this would yield 
9. Applying this over all edges of the tree recursively, we 
can find the exact phase differences between all neighboring 
nodes and since the graph is connected we can find the 
phase of any j £ S hy adding up the phase differences over 
all edges on the path between the initial node and j. This 
can be done for all nodes in 0{k^) time via DFS. Overall, 
the original signal x can be found up to an overall phase 
ambiguity in 0{mn) time. 



